Tarea pronósticos 05

library("easypackages")
packages("tidyverse","fpp3", "tsibble", "feasts","fable", "patchwork")
library("USgas")

Flujo de trabajo limpio de pronóstico

• Preparación de los datos (limpieza)
• Gráfica de los datos (visualización)
• Definición del modelo (especificación)
• Entrenamiento del modelo (estimación)
• Revisar el desempeño del modelo (evaluación)
• Producir pronósticos

Preparación de los datos (limpieza)

Cargar datos en R, o incluso omitir los datos vacíos (NA), filtrado de la serie (mediante paquetería Tsibble).

Gráfica de los datos (visualización)

global_economy %>%
  filter(Country == "Sweden") %>%
  autoplot(GDP) +
    ggtitle("PIB de Suecia") + ylab("$US billions")

Definición del modelo (especificación)

Describir el modelo. Hay distintos tipos, se usará por ejemplo, un modelo lineal de series de tiempo con TLSM:

TSLM(GDP  ~ trend())

<TSLM model definition>

Entrenamiento del modelo (estimación)

Una vez se define el modelo, lo que sigue es entrenar al modelo. Lo que significa pasarle los datos para que, estdísticamente, encuentre los parámetros que realizan el mejor ajuste posible.

fit <- global_economy %>%
  model(Modelo_tendencia = TSLM(GDP ~ trend()))

Warning: 7 errors (1 unique) encountered for Modelo_tendencia
[7] 0 (non-NA) cases

Ahora tenemos un modelo lineal ajustado y el objeto resultante es un mable (model table).

fit

# A mable: 263 x 2
# Key:     Country [263]
   Country             Modelo_tendencia
   <fct>                        <model>
 1 Afghanistan                   <TSLM>
 2 Albania                       <TSLM>
 3 Algeria                       <TSLM>
 4 American Samoa                <TSLM>
 5 Andorra                       <TSLM>
 6 Angola                        <TSLM>
 7 Antigua and Barbuda           <TSLM>
 8 Arab World                    <TSLM>
 9 Argentina                     <TSLM>
10 Armenia                       <TSLM>
# … with 253 more rows

Revisar el desempeño del modelo (evaluación)

Ver como se ajusta el modelo a los datos y compararlo con otros modelos.

Producir pronósticos

Una vez tenemos un modelo ajustado, podemos hacer pronósticos. Usamos forecast(), h = periodo a pronosticar

fcst <- fit %>% forecast(h = 36) # h = "3 years"
fcst

# A fable: 9,468 x 5 [1Y]
# Key:     Country, .model [263]
   Country     .model            Year                 GDP        .mean
   <fct>       <chr>            <dbl>              <dist>        <dbl>
 1 Afghanistan Modelo_tendencia  2018 N(1.6e+10, 1.3e+19) 16205101654.
 2 Afghanistan Modelo_tendencia  2019 N(1.7e+10, 1.3e+19) 16511878141.
 3 Afghanistan Modelo_tendencia  2020 N(1.7e+10, 1.3e+19) 16818654627.
 4 Afghanistan Modelo_tendencia  2021 N(1.7e+10, 1.3e+19) 17125431114.
 5 Afghanistan Modelo_tendencia  2022 N(1.7e+10, 1.3e+19) 17432207600.
 6 Afghanistan Modelo_tendencia  2023 N(1.8e+10, 1.3e+19) 17738984086.
 7 Afghanistan Modelo_tendencia  2024 N(1.8e+10, 1.3e+19) 18045760573.
 8 Afghanistan Modelo_tendencia  2025 N(1.8e+10, 1.3e+19) 18352537059.
 9 Afghanistan Modelo_tendencia  2026 N(1.9e+10, 1.3e+19) 18659313546.
10 Afghanistan Modelo_tendencia  2027 N(1.9e+10, 1.3e+19) 18966090032.
# … with 9,458 more rows

Métodos sencillos de pronóstico

Conocidos como benchmark son básicos pero pueden ser muy utiles.

bricks <- aus_production %>% filter_index("1970" ~ "2004")
bricks

# A tsibble: 137 x 7 [1Q]
   Quarter  Beer Tobacco Bricks Cement Electricity   Gas
     <qtr> <dbl>   <dbl>  <dbl>  <dbl>       <dbl> <dbl>
 1 1970 Q1   387    6807    386   1049       12328    12
 2 1970 Q2   357    7612    428   1134       14493    18
 3 1970 Q3   374    7862    434   1229       15664    23
 4 1970 Q4   466    7126    417   1188       13781    20
 5 1971 Q1   410    7255    385   1058       13299    19
 6 1971 Q2   370    8076    433   1209       15230    23
 7 1971 Q3   379    8405    453   1199       16667    28
 8 1971 Q4   487    7974    436   1253       14484    24
 9 1972 Q1   419    6500    399   1070       13838    24
10 1972 Q2   378    7119    461   1282       15919    34
# … with 127 more rows

bricks %>% autoplot(Bricks)

Método del promedio (media)

En este pronósitco, las predicciones de todos los valores futuros son la media de los datos históricos.

bricks %>% model(MEAN(Bricks))

# A mable: 1 x 1
  `MEAN(Bricks)`
         <model>
1         <MEAN>

Método ingenuo (Naive method)

Se toma el último valor como el pronóstico para todos los valores futuros. También se les conoce como pronósitcos de caminata aleatoria, ya que son óptimos en una serie que siga una.

bricks %>% model(NAIVE(Bricks),
                 RW(Bricks)) # hace exactamente lo mismo que NAIVE()

# A mable: 1 x 2
  `NAIVE(Bricks)` `RW(Bricks)`
          <model>      <model>
1         <NAIVE>      <NAIVE>

Método ingenuo estacional (seasonal Naive)

Se le agrega un componente a Naive para lidiar con datos altamente estacionales.

bricks %>% model(SNAIVE(Bricks))

# A mable: 1 x 1
  `SNAIVE(Bricks)`
           <model>
1         <SNAIVE>

vic_elec %>% 
  autoplot(Demand) %>% 
  plotly::ggplotly()

vic_elec %>% 
  model(SNAIVE(Demand)) %>% 
  forecast(h = "1 day") %>% 
  autoplot(vic_elec %>% filter_index("2014-12-30" ~ .))

Warning: There was 1 warning in `filter()`.
ℹ In argument: `time_in(Time, ...)`.
Caused by warning:
! Argument 'roll' is deprecated. Deprecated in version '1.8.4'.

Método de la deriva (drift)

Variación de Naive que permite que aumente o disminuya el pronótico en el tiempo. El aumento del cambio es el cambio promedio en los datos históricos.

bricks %>% model(
  RW(Bricks ~ drift())
)

# A mable: 1 x 1
  `RW(Bricks ~ drift())`
                 <model>
1          <RW w/ drift>

Ejemplos

# Set training data from 1992 to 2006
train <- aus_production %>% filter_index("1992 Q1" ~ "2006 Q4")
# Fit the models
beer_fit <- train %>%
  model(
    Mean             = MEAN(Beer),
    `Naïve`          = NAIVE(Beer),
    `Seasonal naïve` = SNAIVE(Beer),
    Drift            = RW(Beer ~ drift())
  )
# Generate forecasts for 14 quarters
beer_fc <- beer_fit %>% forecast(h= 14)

# Plot forecasts against actual values
beer_fc %>%
  autoplot(filter_index(aus_production, "1992 Q1" ~ .), level = NULL) +
  ggtitle("Forecasts for quarterly beer production") +
  xlab("Year") + ylab("Megalitres") +
  guides(colour=guide_legend(title="Forecast")) +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick",
             linetype = "dashed") +
  annotate("label", x = c(as.Date("2003-01-01"),as.Date("2009-01-01")),
           y = 550, label = c("Train set", "Test set"),
           color = c("black","blue"))

beer_fc %>%
  filter(.model == "Seasonal naïve") %>% 
  autoplot(filter_index(aus_production, "1992 Q1" ~ .)) +
  ggtitle("Seasonal naive forecast for quarterly beer production") +
  xlab("Quarter") + ylab("Megalitres") +
  geom_vline(xintercept = as.Date("2007-01-01"), color = "firebrick", linetype = "dashed") +
  annotate("label", x = c(as.Date("2003-01-01"), as.Date("2009-01-01")), y = 550, label = c("Train set", "Test set"), color = c("black", "blue"))

Otro ejemplo:

gafa_stock %>% distinct(Symbol)

# A tibble: 4 × 1
  Symbol
  <chr> 
1 AAPL  
2 AMZN  
3 FB    
4 GOOG  

# Re-index based on trading days
google_stock <- gafa_stock %>%
  filter(Symbol == "GOOG") %>%
  mutate(day = row_number()) %>%
  update_tsibble(index = day, regular = TRUE)

# Filter the year of interest
google_2015 <- google_stock %>% filter(year(Date) == 2015)
# Fit the models
google_fit <- google_2015 %>%
  model(
    Mean    = MEAN(Close),
    `Naïve` = NAIVE(Close),
    Drift   = NAIVE(Close ~ drift()),
    SNAIVE  = SNAIVE(Close)
  )

Warning: 1 error encountered for SNAIVE
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.

google_fit

# A mable: 1 x 5
# Key:     Symbol [1]
  Symbol    Mean   Naïve         Drift       SNAIVE
  <chr>  <model> <model>       <model>      <model>
1 GOOG    <MEAN> <NAIVE> <RW w/ drift> <NULL model>

# Produce forecasts for the 19 trading days in January 2015
google_fc <- google_fit %>% forecast(h = 19)

google_fc

# A fable: 76 x 5 [1]
# Key:     Symbol, .model [4]
   Symbol .model   day        Close .mean
   <chr>  <chr>  <dbl>       <dist> <dbl>
 1 GOOG   Mean     505 N(602, 6766)  602.
 2 GOOG   Mean     506 N(602, 6766)  602.
 3 GOOG   Mean     507 N(602, 6766)  602.
 4 GOOG   Mean     508 N(602, 6766)  602.
 5 GOOG   Mean     509 N(602, 6766)  602.
 6 GOOG   Mean     510 N(602, 6766)  602.
 7 GOOG   Mean     511 N(602, 6766)  602.
 8 GOOG   Mean     512 N(602, 6766)  602.
 9 GOOG   Mean     513 N(602, 6766)  602.
10 GOOG   Mean     514 N(602, 6766)  602.
# … with 66 more rows

# A better way using a tsibble to determine the forecast horizons
google_jan_2016 <- google_stock %>%
  filter(yearmonth(Date) == yearmonth("2016 Jan"))
google_fc <- google_fit %>% forecast(google_jan_2016)
# Plot the forecasts
google_fc %>%
  autoplot(google_2015, level = NULL) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))

Warning: Removed 19 rows containing missing values (`()`).

p1 <- google_fc %>%
  filter(.model == "Drift") %>% 
  autoplot(google_2015) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))

p2 <- google_fc %>%
  filter(.model == "Naïve") %>% 
  autoplot(google_2015) +
    autolayer(google_jan_2016, Close, color='black') +
    ggtitle("Google stock (daily ending 31 Dec 2015)") +
    xlab("Day") + ylab("Closing Price (US$)") +
    guides(colour=guide_legend(title="Forecast"))

p1 / p2

Valores ajustados (fitted) y residuales

• Valores reales: Datos históricos de la serie de tiempo.
• Valores ajustados: Valores (fitted) que tratan de pronosticar los valores reales
• Residuales: Valores que el modelo no logró capturar. la diferencia entre valores reales y ajustados.

augment() obtiene los valores ajustados y residuales de un modelo:

augment(beer_fit)

# A tsibble: 240 x 6 [1Q]
# Key:       .model [4]
   .model Quarter  Beer .fitted .resid .innov
   <chr>    <qtr> <dbl>   <dbl>  <dbl>  <dbl>
 1 Mean   1992 Q1   443    436.   6.55   6.55
 2 Mean   1992 Q2   410    436. -26.4  -26.4 
 3 Mean   1992 Q3   420    436. -16.4  -16.4 
 4 Mean   1992 Q4   532    436.  95.6   95.6 
 5 Mean   1993 Q1   433    436.  -3.45  -3.45
 6 Mean   1993 Q2   421    436. -15.4  -15.4 
 7 Mean   1993 Q3   410    436. -26.4  -26.4 
 8 Mean   1993 Q4   512    436.  75.6   75.6 
 9 Mean   1994 Q1   449    436.  12.6   12.6 
10 Mean   1994 Q2   381    436. -55.4  -55.4 
# … with 230 more rows

Diagnóstico de residuales

Es importante analizar los residuales de un modelo. Si los residuales presentan patrones, es signo de que el modelo se puede mejorar.

Los residuales de un buen modelo deben de:

• No estar autocorrelacionadas. Si se detectan correlaciones entre residuos,todavía hay información útil que se debe modelar.
• La media de los residuos es cero. Si es distinta a cero, el pronóstico está sesgado.

Opcionales: - Los residuos tienen una varianza constante. - Los residuos se distribuyen de manera normal.

Transformaciones Box-Cox pueden ayudar a mejorar los residuales

google_2015 %>% autoplot(Close) +
  xlab("Day") + ylab("Closing Price (US$)") +
  ggtitle("Google Stock in 2015")

aug <- google_2015 %>% 
  model(NAIVE(Close)) %>% 
  augment()

aug

# A tsibble: 252 x 7 [1]
# Key:       Symbol, .model [1]
   Symbol .model         day Close .fitted  .resid  .innov
   <chr>  <chr>        <int> <dbl>   <dbl>   <dbl>   <dbl>
 1 GOOG   NAIVE(Close)   253  522.     NA   NA      NA    
 2 GOOG   NAIVE(Close)   254  511.    522. -10.9   -10.9  
 3 GOOG   NAIVE(Close)   255  499.    511. -11.8   -11.8  
 4 GOOG   NAIVE(Close)   256  498.    499.  -0.855  -0.855
 5 GOOG   NAIVE(Close)   257  500.    498.   1.57    1.57 
 6 GOOG   NAIVE(Close)   258  493.    500.  -6.47   -6.47 
 7 GOOG   NAIVE(Close)   259  490.    493.  -3.60   -3.60 
 8 GOOG   NAIVE(Close)   260  493.    490.   3.61    3.61 
 9 GOOG   NAIVE(Close)   261  498.    493.   4.66    4.66 
10 GOOG   NAIVE(Close)   262  499.    498.   0.915   0.915
# … with 242 more rows

# aug %>% 
#   features(.resid, mean(.,na.rm = TRUE))

aug %>% pull(.resid) %>% mean(na.rm = TRUE) 

[1] 0.9439931

aug %>% autoplot(.resid) + xlab("Día") + ylab("") +
  ggtitle("Residuales del método naïve")

Warning: Removed 1 row containing missing values (`geom_line()`).

Parece que la media es casi cero y la variación es constante, excepto un outlier.

aug %>%
  ggplot(aes(x = .resid)) +
  geom_histogram() +
  ggtitle("Histograma de los residuales")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

aug %>% ACF(.resid)

# A tsibble: 23 x 4 [1]
# Key:       Symbol, .model [1]
   Symbol .model            lag      acf
   <chr>  <chr>        <cf_lag>    <dbl>
 1 GOOG   NAIVE(Close)        1  0.0976 
 2 GOOG   NAIVE(Close)        2 -0.0726 
 3 GOOG   NAIVE(Close)        3 -0.0748 
 4 GOOG   NAIVE(Close)        4 -0.0433 
 5 GOOG   NAIVE(Close)        5 -0.0398 
 6 GOOG   NAIVE(Close)        6  0.0206 
 7 GOOG   NAIVE(Close)        7 -0.0652 
 8 GOOG   NAIVE(Close)        8 -0.0291 
 9 GOOG   NAIVE(Close)        9 -0.0379 
10 GOOG   NAIVE(Close)       10 -0.00687
# … with 13 more rows

aug %>% ACF(.resid) %>% autoplot() + ggtitle("ACF of residuals")

Se puede ver que los residuos no están autocorrelacionados (ACF), por lo tanto este método (NAIVE) es bueno para esta serie de tiempo, y los pronósticos derivados de este método pueden ser muy buenos.

Todo en una sola gráfica:

google_2015 %>% 
  model(NAIVE(Close)) %>% 
  gg_tsresiduals()

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

aus_production %>% 
  filter_index("1992" ~ .) %>%
  model(SNAIVE(Beer)) %>% 
  gg_tsresiduals()

Warning: Removed 4 rows containing missing values (`geom_line()`).

Warning: Removed 4 rows containing missing values (`geom_point()`).

Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

aus_production %>% 
  filter_index("1992" ~ .) %>%
  gg_tsdisplay(Beer)

Veamos con el método de Media:

google_2015 %>% 
  model(MEAN(Close)) %>% 
  augment() %>% 
  pull(.resid) %>% 
  mean(na.rm = TRUE)

[1] 2.710222e-14

google_2015 %>% 
  model(MEAN(Close)) %>% 
  gg_tsresiduals() + 
  ggtitle("Diagnóstico de residuales para el modelo de Media")

Como podemos ver, el método de la media no es óptimo para esta serie de acciones, siendo que la media es 0, al graficar los residuales vemos que muestra un patón (igual al de los datos menos la media), existen claramente autocorrelaciones de una caminata aleatoria y no hay ninguna distribución normal. estos tres factores lo hacen un método nada efectivo para acciones.

Tests de Portmanteau de autocorrelación

Tests para analizar de una forma más formal la presencia o ausencia de autocorrelación en los residuales. Nos ayuda a ver si las primeras \(h\) autocorrelaciones son significativamente distintas de cero o no.

Test de Box-Pierce

Sumatoria de \(r^2_k\) donde \(h\) es el rezago máximo a considerar y \(T\) es la cantidad de observaciones en la muestra. Si \(r^2_k\) es pequeña, entonces el resultado será pequeño. Se recomienda usar \(h = 10\) para datos no estacionales, \(h = 2m\) para datos estaionales (m es el periodo estacional) o como máximo \(h = T/5\).

Hay que tomar en cuenta, que la prueba no es tan buena cuando \(h\) es grande (mayor a \(T/5\)).

Test de Ljung-Box

Generalmente más preciso que Box-Pierce, de igual forma, valores grandes en el resultado indica que no es ruido blanco y los residuales sí están autocorrelacionados.

Hipótesis nula

La hipótesis nula dice que la serie en cuestión no está autocorrelacionada, o sea, es riudo blanco. Para unos residuales, que la hipótesis nula sea cierta es bueno, ya que significa que no están correlacionados los residuales.

¿Cómo saber si es cierta o no?

Siendo \(\alpha\) el nivel de significancia o nivel máximo de error dispuestos a aceptar (usaremos 5%) y p-value un avlor resultante de las pruebas:

Si p-value es menor a \(\alpha\), entonces rechazamos la hiótesis nula, de lo contrario, la aceptamos.

O sea:

• \(p-value > \alpha\), H0 es cierta, residuales no están autocorrelacionados.
• \(p-value < \alpha\), H0 es falsa, residuales están autocorrelacionados.

Ejemplos:

# lag=h and fitdf=K
aug %>% features(.resid, box_pierce, lag=10, dof=0)

# A tibble: 1 × 4
  Symbol .model       bp_stat bp_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.74     0.654

aug %>% features(.resid, ljung_box, lag=10, dof=0)

# A tibble: 1 × 4
  Symbol .model       lb_stat lb_pvalue
  <chr>  <chr>          <dbl>     <dbl>
1 GOOG   NAIVE(Close)    7.91     0.637

Como podemos ver para el método Naive, el p-value es mayor al 5%, por lo tanto sus residuos no están autocrrelacionados, son ruido blanco.

google_2015 %>% 
  model(MEAN(Close)) %>% augment() %>% 
  features(.resid, box_pierce, lag = 10, dof = 0)

# A tibble: 1 × 4
  Symbol .model      bp_stat bp_pvalue
  <chr>  <chr>         <dbl>     <dbl>
1 GOOG   MEAN(Close)   2024.         0

google_2015 %>% 
  model(MEAN(Close)) %>% augment() %>% 
  features(.resid, ljung_box, lag = 10, dof = 0)

# A tibble: 1 × 4
  Symbol .model      lb_stat lb_pvalue
  <chr>  <chr>         <dbl>     <dbl>
1 GOOG   MEAN(Close)   2083.         0

Para el caso de Google, no se logra distinguir los residuales del ruido blanco.

Intervalos de predicción

\(c\) es el porcentaje de cobertura de probabilidad. Normalmente se usará 80% y 95%.

Al tener incertidumbre en el ponóstico, entre maoyr sea el intervalo, menos preciso será nuestro pronóstico.

Intervalos de predicción de un paso

Cunado ser realizan pronósitocs de un paso, la desviación estándar del pronóstico es prácticamente la misma que la desviación estándar de los residuos.

Intervalos de predicción de paso múltiple

Entre más aumenta el horizonte de pronóstico, mayor será el intervalo de pronóstico, Esto es, \(\simga_h\) incrementa con \(h\), entonces necesitamos estimaciones de \(\sigma_h\).

Métodos de referencia

Usando fable:

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  hilo()

# A tsibble: 10 x 7 [1]
# Key:       Symbol, .model [1]
   Symbol .model         day        Close .mean                  `80%`
   <chr>  <chr>        <dbl>       <dist> <dbl>                 <hilo>
 1 GOOG   NAIVE(Close)   505  N(759, 125)  759. [744.5400, 773.2200]80
 2 GOOG   NAIVE(Close)   506  N(759, 250)  759. [738.6001, 779.1599]80
 3 GOOG   NAIVE(Close)   507  N(759, 376)  759. [734.0423, 783.7177]80
 4 GOOG   NAIVE(Close)   508  N(759, 501)  759. [730.1999, 787.5601]80
 5 GOOG   NAIVE(Close)   509  N(759, 626)  759. [726.8147, 790.9453]80
 6 GOOG   NAIVE(Close)   510  N(759, 751)  759. [723.7543, 794.0058]80
 7 GOOG   NAIVE(Close)   511  N(759, 876)  759. [720.9399, 796.8202]80
 8 GOOG   NAIVE(Close)   512 N(759, 1002)  759. [718.3203, 799.4397]80
 9 GOOG   NAIVE(Close)   513 N(759, 1127)  759. [715.8599, 801.9001]80
10 GOOG   NAIVE(Close)   514 N(759, 1252)  759. [713.5329, 804.2272]80
# … with 1 more variable: `95%` <hilo>

google_2015 %>%
  model(NAIVE(Close)) %>%
  forecast(h = 10) %>%
  autoplot(google_2015)

Intervalos de predicción con residuales bootstrap

Cuando no es razonable asumir normlidad en los residuos, podemos aplicarles bootstraping, ya que solo asume la no autocorrelación.

Se obteienen muchos escenarios futuros posibles. Para ver algunos, se usa generate:

fit <- google_2015 %>%
  model(NAIVE(Close))

sim <- fit %>%  generate(h = 30, times = 5, bootstrap = TRUE, seed = 123) # only works with dev version

sim

# A tsibble: 150 x 5 [1]
# Key:       Symbol, .model, .rep [5]
   Symbol .model         day .rep   .sim
   <chr>  <chr>        <dbl> <chr> <dbl>
 1 GOOG   NAIVE(Close)   505 1      744.
 2 GOOG   NAIVE(Close)   506 1      747.
 3 GOOG   NAIVE(Close)   507 1      733.
 4 GOOG   NAIVE(Close)   508 1      737.
 5 GOOG   NAIVE(Close)   509 1      739.
 6 GOOG   NAIVE(Close)   510 1      733.
 7 GOOG   NAIVE(Close)   511 1      728.
 8 GOOG   NAIVE(Close)   512 1      730.
 9 GOOG   NAIVE(Close)   513 1      721.
10 GOOG   NAIVE(Close)   514 1      713.
# … with 140 more rows

Generamos 5 escenarios futuros posibles para los siguientes 30 días de trading.

google_2015 %>%
  ggplot(aes(x = day)) +
  geom_line(aes(y = Close), size = 1) +
  geom_line(aes(y = .sim, colour = as.factor(.rep)), data = sim, size = 1) +
  ggtitle("Google closing stock price") +
  guides(col = FALSE)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.

Con esto, podemos obtener intervalos de predicción, al calcular los percentiles de los escenarios futuros. El resultado es el ntervalo de predicción bootstrapped. Esto se puede lograr con forecast:

fc <- fit %>% forecast(h = 30, bootstrap = TRUE)
fc

# A fable: 30 x 5 [1]
# Key:     Symbol, .model [1]
   Symbol .model         day        Close .mean
   <chr>  <chr>        <dbl>       <dist> <dbl>
 1 GOOG   NAIVE(Close)   505 sample[5000]  759.
 2 GOOG   NAIVE(Close)   506 sample[5000]  759.
 3 GOOG   NAIVE(Close)   507 sample[5000]  759.
 4 GOOG   NAIVE(Close)   508 sample[5000]  759.
 5 GOOG   NAIVE(Close)   509 sample[5000]  759.
 6 GOOG   NAIVE(Close)   510 sample[5000]  759.
 7 GOOG   NAIVE(Close)   511 sample[5000]  759.
 8 GOOG   NAIVE(Close)   512 sample[5000]  759.
 9 GOOG   NAIVE(Close)   513 sample[5000]  759.
10 GOOG   NAIVE(Close)   514 sample[5000]  759.
# … with 20 more rows

fc %>% autoplot(google_2015) +
  ggtitle("Google closing stock price")

fc <- fit %>% forecast(h = 30)
fc

# A fable: 30 x 5 [1]
# Key:     Symbol, .model [1]
   Symbol .model         day        Close .mean
   <chr>  <chr>        <dbl>       <dist> <dbl>
 1 GOOG   NAIVE(Close)   505  N(759, 125)  759.
 2 GOOG   NAIVE(Close)   506  N(759, 250)  759.
 3 GOOG   NAIVE(Close)   507  N(759, 376)  759.
 4 GOOG   NAIVE(Close)   508  N(759, 501)  759.
 5 GOOG   NAIVE(Close)   509  N(759, 626)  759.
 6 GOOG   NAIVE(Close)   510  N(759, 751)  759.
 7 GOOG   NAIVE(Close)   511  N(759, 876)  759.
 8 GOOG   NAIVE(Close)   512 N(759, 1002)  759.
 9 GOOG   NAIVE(Close)   513 N(759, 1127)  759.
10 GOOG   NAIVE(Close)   514 N(759, 1252)  759.
# … with 20 more rows

fc %>% autoplot(google_2015) +
  ggtitle("Google closing stock price")

Pronósticos con transformaciones

Si se le hace un pronóstico a una serie anteriormente transformada (Box-cox por ejemplo), al terminar el pronóstico, se tiene que hacer una transformación inversa para ver la versión original del pronóstico.

Transformación inversa de Box-Cox:

\(y_t=\)

\({exp(w_t), \text{ si }\lambda = 0}.\)

\((\lambda w_t^\lambda+1)^{1/\lambda}, \text{ en otro caso}.\)

La paquetería fable realiza la transofrmación ivnersa en automático, cuando se especifica en la definición del modelo.

Ajustes por sesgo

Un problema al realizar transformaciones matemáticas, como Box–Cox es que las estimaciones puntuales re transformadas ya no presentan la media de la distribución de predicción, sino que la mediana. Esto puede llegar a ser problema dependiendo del contexto, pero normalmente no causa mucha diferencia.

Ajustar por sesgo es la transformación inversa de la media, entre más grande sea la varianza, mayor será la diferencia entre la media y la mediana.

eggs <- as_tsibble(fma::eggs)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

eggs %>% 
  model(RW(log(value) ~ drift())) %>% 
  forecast(h=50) %>% 
  autoplot(eggs, level = 80,  
           point_forecast = lst(mean, median))

Warning in ggplot2::geom_point(mapping = mapping, data =
dplyr::semi_join(object, : Ignoring unknown aesthetics: linetype

Por lo tanto, para tener un pronóstico sesgado, se usa point_forecast = lst(median)

Pronósticos con descomposición

La descomposición de series de tiempo puede ser útil para producir pronósticos.

El pronóstico se realiza en dos pasos: - Un pronóstico para el componente estacional. - Un pronósitco separado para la serie desestacionalizada.

En realidad, SNAIVE es el pronóstido de solo el componente estacional.

(us_retail_employment <- fpp3::us_employment %>%
  filter(year(Month) >= 1990, Title == "Retail Trade"))

# A tsibble: 357 x 4 [1M]
# Key:       Series_ID [1]
      Month Series_ID     Title        Employed
      <mth> <chr>         <chr>           <dbl>
 1 1990 Jan CEU4200000001 Retail Trade   13256.
 2 1990 Feb CEU4200000001 Retail Trade   12966.
 3 1990 Mar CEU4200000001 Retail Trade   12938.
 4 1990 Apr CEU4200000001 Retail Trade   13012.
 5 1990 May CEU4200000001 Retail Trade   13108.
 6 1990 Jun CEU4200000001 Retail Trade   13183.
 7 1990 Jul CEU4200000001 Retail Trade   13170.
 8 1990 Aug CEU4200000001 Retail Trade   13160.
 9 1990 Sep CEU4200000001 Retail Trade   13113.
10 1990 Oct CEU4200000001 Retail Trade   13185.
# … with 347 more rows

us_retail_employment %>% 
autoplot(Employed)

dcmp <- us_retail_employment %>%
  model(STL(Employed ~ trend(window = 7), robust=TRUE)) %>%
  components() %>%
  select(-.model)

dcmp

# A tsibble: 357 x 7 [1M]
# Key:       Series_ID [1]
   Series_ID        Month Employed  trend season_year remainder season_adjust
   <chr>            <mth>    <dbl>  <dbl>       <dbl>     <dbl>         <dbl>
 1 CEU4200000001 1990 Jan   13256. 13326.       -72.6     2.85         13328.
 2 CEU4200000001 1990 Feb   12966. 13297.      -273.    -58.1          13239.
 3 CEU4200000001 1990 Mar   12938. 13270.      -295.    -36.3          13233.
 4 CEU4200000001 1990 Apr   13012. 13230.      -216.     -1.59         13228.
 5 CEU4200000001 1990 May   13108. 13223.      -119.      4.72         13228.
 6 CEU4200000001 1990 Jun   13183. 13212.       -30.6     0.978        13213.
 7 CEU4200000001 1990 Jul   13170. 13198.       -32.8     5.23         13203.
 8 CEU4200000001 1990 Aug   13160. 13178.       -20.1     1.84         13180.
 9 CEU4200000001 1990 Sep   13113. 13149.       -45.5     9.76         13159.
10 CEU4200000001 1990 Oct   13185. 13113.        64.2     8.28         13121.
# … with 347 more rows

dcmp %>%
  model(NAIVE(season_adjust)) %>%
  forecast() %>%
  autoplot(dcmp) + ylab("New orders index") +
  ggtitle("Pronóstico naïve de los datos desestacionalizados")

Le podemos agregar nuevamente la estacionalidad con la función decomposition_model():

us_retail_employment %>%
  model(stlf = decomposition_model(
             STL(Employed ~ trend(window = 7), robust = TRUE),
             NAIVE(season_adjust)
  )) %>%
  forecast() %>%
  autoplot(us_retail_employment)

Evaluación del desempeño de los pronósticos

Conjuntos de entrenamiento y prueba

El tamaño de la prueba es, generalmente, del 20% del total de datos disponibles.

• Un modelo que se ajusta muy bien a los datos de entrenamiento no necesariamente produce los mejore spronósticos.
• Al aumentar los parámetros, podemos llegar a un ajuste perfecto del modelo a los datos.
• Puede darse un efecto de sobre ajuste (over-fitting) y esto es tan malo como tener un muy mal ajuste.

Funciones para sementar las series de tiempo

top_n nos permite obtener las n observaciones más extremas.

gafa_stock %>%
  group_by(Symbol) %>%
  top_n(1, Close)

# A tsibble: 4 x 8 [!]
# Key:       Symbol [4]
# Groups:    Symbol [4]
  Symbol Date        Open  High   Low Close Adj_Close   Volume
  <chr>  <date>     <dbl> <dbl> <dbl> <dbl>     <dbl>    <dbl>
1 AAPL   2018-10-03  230.  233.  230.  232.      230. 28654800
2 AMZN   2018-09-04 2026. 2050. 2013  2040.     2040.  5721100
3 FB     2018-07-25  216.  219.  214.  218.      218. 58954200
4 GOOG   2018-07-26 1251  1270. 1249. 1268.     1268.  2405600

Errores de pronóstico

Diferencia entre el valor real ocurrido y el dato pronosticado.

Los errores de pronóstico son distintos de los residuales en:

• Los residuales se calculan con los datos de entrenamiento, mientras que los errores de pronóstico se calculan con los de prueba.
• Los residuos se calculan mediante pronósticos de un paso, los errores pueden ser multi-step.

Errores dependientes de la escala de los datos

Dependen de la escala de la serie de tiempo, por lo que no se puede comparar con series de tiempo con otras unidades.

Los dos más utilizados son:

• MAE, Mean absolute error, el promedio del error.
• RMSE, Root mean squared error, la raíz cuadrada del promedio del error.

MAE es muy utilizado por su facilidad de cómputo e interpretación.

Errores porcentuales

Son un porcentaje y son muy utilizados para comparar el desempeño de pronósticos de distintos conjuntos de datos.

-MAPE, Mean absolute percentage error, el promedio del error.

El problema es que si un dato es 0, este error será infinito, por lo tanto se inventó:

-sMAPE, symmetric MAPE, resuelve el problema del MAPE.

Sin embargo, no es recomendable en la práctica.

Errores escalados

ALternativa para porcentuales para comparar con otras series.

• MASE, Mean absolute scalated error, promedio de los errores escalados.

recent_production <- aus_production %>% filter(year(Quarter) >= 1992)
beer_train <- recent_production %>% filter(year(Quarter) <= 2007)

beer_fit <- beer_train %>%
  model(
    Mean = MEAN(Beer),
    `Naïve` = NAIVE(Beer),
    `Seasonal naïve` = SNAIVE(Beer),
    Drift = RW(Beer ~ drift())
  )

beer_fc <- beer_fit %>%
  forecast(h = 10)

beer_fc %>%
  autoplot(recent_production, level = NULL) +
  xlab("Year") + ylab("Megalitres") +
  ggtitle("Forecasts for quarterly beer production") +
  guides(colour=guide_legend(title="Forecast"))

Usamos accuracy para ver los errores en el ajuste de modelos (mabble) a los datos de entrenamiento:

accuracy(beer_fit)

# A tibble: 4 × 10
  .model         .type           ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
  <chr>          <chr>        <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
1 Mean           Training  0         43.6  35.2 -0.937  7.89  2.46  2.60 -0.109
2 Naïve          Training  4.76e- 1  65.3  54.7 -0.916 12.2   3.83  3.89 -0.241
3 Seasonal naïve Training -2.13e+ 0  16.8  14.3 -0.554  3.31  1     1    -0.288
4 Drift          Training -5.41e-15  65.3  54.8 -1.03  12.2   3.83  3.89 -0.241

Ahora la usamos para ver los errores de pronóstico (fabble) comparado a los datos de prueba:

beer_fc |>  accuracy(recent_production)

# A tibble: 4 × 10
  .model         .type    ME  RMSE   MAE    MPE  MAPE  MASE RMSSE    ACF1
  <chr>          <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl>
1 Drift          Test  -54.0  64.9  58.9 -13.6  14.6  4.12  3.87  -0.0741
2 Mean           Test  -13.8  38.4  34.8  -3.97  8.28 2.44  2.29  -0.0691
3 Naïve          Test  -51.4  62.7  57.4 -13.0  14.2  4.01  3.74  -0.0691
4 Seasonal naïve Test    5.2  14.3  13.4   1.15  3.17 0.937 0.853  0.132 

Los errores de ajuste (datos de entrenamiento) sirven para ver que modelos utilizar para pronosticar.

Los errores de pronóstico (datos de prueba) nos da mayor claridad de cuáles modelos parecen producir los mejores pronósticos.

Después de esto se recalculan los modelos utilizando toda la información (sin separar en conjuntos de entrenamiento y prueba) y se producen los pronósticos reales a presentar.

Tarea

Series:

1. Serie de ventas de ropa en el Estado de Victoria, Australia

(clothes_victoria <- filter(aus_retail, State == "Victoria", Industry == "Clothing retailing") |> select(Month, Turnover))

# A tsibble: 441 x 2 [1M]
      Month Turnover
      <mth>    <dbl>
 1 1982 Apr     93.6
 2 1982 May     95.3
 3 1982 Jun     85.2
 4 1982 Jul     91.6
 5 1982 Aug     85.2
 6 1982 Sep     89.5
 7 1982 Oct     93  
 8 1982 Nov    108. 
 9 1982 Dec    148. 
10 1983 Jan     81.6
# … with 431 more rows

train <- filter_index(clothes_victoria, "2000 Jan" ~ "2015 Dec")
autoplot(clothes_victoria, Turnover)

fit <- model(train, 
             Mean = MEAN(Turnover), 
             Naive = NAIVE(Turnover), 
             Snaive = SNAIVE(Turnover), 
             Drift = RW(Turnover ~ drift()))

fc <- forecast(fit, h = "3 years")

fc |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .), level = NULL) +
  ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
  xlab("Mes") + ylab("MegaWhats/hora") + 
  guides(color=guide_legend(title="Tipo")) +
  geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")

select(fit, Mean) |> gg_tsresiduals() + ggtitle("Media") 

La media no es buena idea, no muestra distribución normal del todo, están autocorrelacionados y la gráfica muestra un patrón.

select(fit, Naive) |> gg_tsresiduals() + ggtitle("Naive") 

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Naive parece ser la mejor opción, tiene distribución normal, su varianza es algo constante y no parece estar muy correlacionado.

select(fit, Snaive) |> gg_tsresiduals() + ggtitle("Snaive") 

Warning: Removed 12 rows containing missing values (`geom_line()`).

Warning: Removed 12 rows containing missing values (`geom_point()`).

Warning: Removed 12 rows containing non-finite values (`stat_bin()`).

Snaive muestra distribución normal, sin embargo, parece presentar cierta correlación y patrón. A pesar de esto, Snaive sigue siendo un buen modelo al ver el gráfico.

select(fit, Drift) |> gg_tsresiduals() + ggtitle("Drift") 

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Drift tiene residuales muy similares a Naive, pueden ser buenas bases para aplicar otros modelos, sin embargo, tomando solo en cuenta los benchmark, probablemente Snaive será el que se va a utilizar.

fc <- filter(fc, .model == "Snaive")
fc_bs <- forecast(fit, h ="3 years", bootstrap = TRUE)
fc_bs <- filter(fc_bs, .model == "Snaive")

fc |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .)) +
  ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
  xlab("Mes") + ylab("MegaWhats/hora") + 
  guides(color=guide_legend(title="Tipo")) +
  geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")

fc_bs |> autoplot(filter_index(clothes_victoria, "2000 Jan" ~ .)) +
  ggtitle("Pronóstico de ventas mensuales de ropa en Victoria") +
  xlab("Mes") + ylab("MegaWhats/hora") + 
  guides(color=guide_legend(title="Tipo")) +
  geom_vline(xintercept = as.Date("2016-01-01"), color = "firebrick", linetype = "dashed")

Al no ser datos muy aleatorios o los datos de una acción, y que los horizontes de pronóstico no son muy diferentes, no es necesario usar Bootstrap.

2. Serie de PIB per capita de Australia

(aus_gdp <- filter(global_economy, Code == "AUS") %>% 
   mutate(GDP_pc = GDP / Population) |> select(Year, GDP_pc))

# A tsibble: 58 x 2 [1Y]
    Year GDP_pc
   <dbl>  <dbl>
 1  1960  1807.
 2  1961  1874.
 3  1962  1851.
 4  1963  1964.
 5  1964  2128.
 6  1965  2277.
 7  1966  2340.
 8  1967  2576.
 9  1968  2719.
10  1969  2986.
# … with 48 more rows

train <- filter_index(aus_gdp, "1970" ~ "2010")
autoplot(aus_gdp, GDP_pc)

fit <- model(train, 
             Mean = MEAN(GDP_pc), 
             Naive = NAIVE(GDP_pc), 
             Snaive = SNAIVE(GDP_pc), 
             Drift = RW(GDP_pc ~ drift()))

Warning: 1 error encountered for Snaive
[1] Non-seasonal model specification provided, use RW() or provide a different lag specification.

fc <- forecast(fit, h = "7 years")

fc |> autoplot(filter_index(aus_gdp, "1970" ~ .), level = NULL) +
  ggtitle("Pronóstico de PIB per Capita de Australia") +
  xlab("Año") + ylab("USD") + 
  guides(color=guide_legend(title="Tipo")) +
  geom_vline(xintercept = 2011, color = "firebrick", linetype = "dashed")

Warning: Removed 7 rows containing missing values (`()`).

Como podemos ver, Snaive no funciona al no ser estacional (se usan datos anuales)

select(fit, Mean) |> gg_tsresiduals() + ggtitle("Media") 

La media no es buena idea, no muestra distribución normal, están autocorrelacionados y la gráfica muestra un patrón muy similar a los datos.

select(fit, Naive) |> gg_tsresiduals() + ggtitle("Naive") 

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Naive parece ser la mejor opción, tiene distribución normal, y a pesar de no tener una varianza muy constante, es la mejor, y no parece estar muy correlacionado.

select(fit, Drift) |> gg_tsresiduals() + ggtitle("Drift") 

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

Drift tiene residuales muy similares a Naive, y tomando en cuenta las gráficas, se decide usar Drift para capturar la tendencia alcista que tiene.

Al no ser datos muy aleatorios o los datos de una acción, y que los horizontes de pronóstico no son muy diferentes, no es necesario usar Bootstrap.